      program bayesgrd5
      implicit real*8 (a-h,l,o-z)
      parameter(nplay=101,ndec=12,ngrid=51)
c
      character*24 fdate
      dimension bgrid(ngrid),pz(ngrid),w0(ngrid),w(ngrid,ngrid)
      dimension ti(nplay,ngrid,ngrid),tl(ngrid,ngrid),pl(ngrid,ngrid)
      dimension fi(nplay,ngrid,ngrid)
      dimension px(ndec,2),pxm(ngrid,ngrid,ndec,2)
      dimension is(nplay,ndec)
      common /mats/is
      common /pmatrx/pxm
c
      open(unit=11,file='ra2011.data')   !with card sleeves on boards
c      open(unit=20,file='bayesgrd5.dat')
c
      read(11,*)((is(i,n),n=1,ndec),i=1,nplay)
c
c-----------------------------------------------------------------
      write(*,'('' ***** bayesgrd5 ******'',20x,a24,
     &//"    Risk & Ambiquity Experiment")')fdate()
      write(*,'(/" All 2011 Sessions"//)')
c
      do 5 i=1,ngrid
         if ((i.eq.1).or.(i.eq.ngrid)) then
            w0(i) = 0.5d0
         else
            w0(i) = 1.d0
         endif
    5 continue
c
      dz=1.d0/float(ngrid-1)
      do 20 i=1,ngrid
         beta = dz*(i-1)
         bgrid(i) = beta
         do 15 j=1,ngrid
            if (j.lt.ngrid) then
               pz(j) = 0.5d0 + 0.5d0*dz*(j-1)
            else
               pz(j) = 0.999d0
            endif
            gam = 2.d0*dlog(pz(j)/(1.d0-pz(j)))
            call fpx(gam,beta,px)
            do 10 k=1,ndec
               pxm(i,j,k,1) = px(k,1)
               pxm(i,j,k,2) = px(k,2)
   10       continue
            pl(i,j) = 0.d0
            tl(i,j) = 0.d0
            w(i,j) = w0(i)*w0(j)
   15    continue
   20 continue
c
      do 50 id=1,nplay
         cumi = 0.d0
         do 35 i=1,ngrid
            do 30 j=1,ngrid
               call fam(id,i,j,t)
               ti(id,i,j) = t
               fi(id,i,j) = t
               cumi = cumi + t*w(i,j)
   30       continue
   35    continue
         do 45 i=1,ngrid
            do 40 j=1,ngrid
               ti(id,i,j) = ti(id,i,j)/cumi
               tl(i,j) = tl(i,j) + ti(id,i,j)
   40       continue
   45    continue
   50 continue
c
      write(*,'(7x,51f7.3)')(bgrid(j),j=1,ngrid)
c      write(20,'(7x,51f7.3)')(bgrid(j),j=1,ngrid)
      xn = 1.d2/float(nplay)
      do 80 id=1,nplay
         cumi = 0.d0
         do 65 i=1,ngrid
            do 60 j=1,ngrid
               pij = tl(i,j) - ti(id,i,j)
               cumi = cumi + w(i,j)*ti(id,i,j)*pij
   60       continue
   65    continue
         do 75 i=1,ngrid
            do 70 j=1,ngrid
               pij = tl(i,j) - ti(id,i,j)
               pij = ti(id,i,j)*pij
               pl(i,j) = pl(i,j) + xn*pij/cumi
   70       continue
   75    continue
   80 continue
c
      xll = 0.d0
      do 100 id=1,nplay
         pid = 0.d0
         do 95 j=1,ngrid
            do 90 k=1,ngrid
               pid = pid + 1.d-2*w(j,k)*pl(j,k)*fi(id,j,k)
   90       continue
   95    continue
         xll = xll + dlog(pid)
  100 continue
      write(*,'(//"LL given g*: ",g14.6/)')xll
c
      gt = 0.d0
      abeta=0.d0
      az=0.d0
      do 120 j0=1,ngrid
         j=ngrid+1-j0
         write(*,'(52f7.3)')pz(j),(pl(i,j),i=1,ngrid)
c         write(20,'(52f7.3)')pz(j),(pl(i,j),i=1,ngrid)
         do 110 i=1,ngrid
            gt = gt + w(i,j)*pl(i,j)
            abeta = abeta + w(i,j)*pl(i,j)*bgrid(i)
            ap = ap + w(i,j)*pl(i,j)*pz(j)
  110    continue
  120 continue
      abeta = abeta/gt
      ap = ap/gt
      agam = 2.d0*dlog(ap/(1.d0-ap))
c
      write(*,'(/"gt = ",f7.3)')gt
      write(*,'(/"abeta, ap, agam: ",3g13.5)')abeta,ap,agam
      end
c***********************************************************************
      subroutine fam(id,i,j,t)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2)
      dimension is(nplay,ndec)
      common /mats/is
      common /pmatrx/pxm
c
      t=1.d0
      do 10 k=1,ndec
         m = is(id,k)
         t = t*pxm(i,j,k,m)
   10 continue
c
      return
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end

